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When twisting a strip of paper or acetate under high longitudinal tension, one 

observes, at some critical load, a buckling of the strip into a regular triangular 

pattern. Very similar triangular facets have recently been observed in solutions to 

QQ , a new set of geometrically-exact equations describing the equilibrium shape of thin 

inextensible elastic strips. Here we formulate a modified boundary-value problem for 

these equations and construct post-buckling solutions in good agreement with the 

fH ' observed pattern in twisted strips. We also study the force-extension and nioment- 

^ |. twist behaviour of these strips by varying the mode number n of triangular facets. 

jy^ ^ Keywords: twisted inextensible strip, triangular buckling pattern, calculus of 

C^ , variations, boundary-value problem, load-displacement behaviour 

CZ2 ' 1. Introduction 

O 

C/3 ■ When twisting a strip of paper or acetate under high longitudinal tension, one 

^~»| observes, at a critical load, a buckling of the strip into a regular triangular pattern 

(see figure [TJ a)). The deformation is reversible. Sheets of paper or acetate are for 

practical purposes inextensible and the observed pattern, consisting of helically 

stacked nearly-flat triangular facets, appears to be nature's way of achieving global 

twisting by means of local bending and minimal stretch. This mode of buckling 

does not appear to have been reported in the literature. Perhaps this is because 

p^ I the phenomenon occurs at relatively large twisting angles and under relatively high 

(T^ . tension (in order to suppress the more common looping instability). The buckling 

C^ ' pattern observed has ridges running at roughly 45° angles to the centreline of the 

1/-^ . strip. The ridges radiate out from vertices on the edge of the strip where stress 

f^ ' concentration occurs. There is superficial similarity with the well-known Yoshimura 

^D . or diamond buckling pattern of thin cylindrical shells (Yoshimura 1930) but this 

pattern requires compressive rather than tensile loading. 

The phenomenon of buckling of an extensible elastic strip under tension and 

twisting moment has been studied before. In the 1930s Green considered buckling 

of twisted strips under constant tensile force, treating both the case of zero and 

C^ ' non-zero end force (Green 1936, Green 1937). Buckling of a twisted orthotropic 

plate into a sinusoidal buckling pattern in the longitudinal direction was studied in 

Crispino & Benson (1986). A numerical investigation of wrinkling of twisted plates 
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was carried in Mockensturm (2001), considering both the case of constant end-to- 
end distance and the case of constant end force, and finding different results in 
the two cases. In the latter case it was found that buckhng may occur in both the 
lateral and longitudinal directions. 

A perturbation method was used recently to further explore the wrinkhng in- 
stabihty under small twist (Coman & Bassom 2008). As the twist increases from 
zero, the surface of a strip first becomes helicoidal, which causes extensional forces 
near the edges while the core domain is in compression. When a critical twist and 
tension are reached, the core domain of the strip buckles into an oscillatory pattern. 
The surface near the edges remains helicoidal; no strain localisation occurs. 

The problem of bending and twisting of an inextensible flat plate has also been 
studied before. Such a plate deforms isometrically and its surface is therefore de- 
velopable (i.e., has zero Gaussian curvature), making the analysis more geomet- 
rical. Sadowsky developed a large-deformation theory of narrow elastic strips as 
early as in 1930 (Sadowsky 1931). Approximate equations for wide strips were de- 
rived in the mid-1950s by Mansfield (see Mansfield 1989). These equations predict 
the distribution of generators of the developable surface while ignoring the actual 
three-dimensional geometry. This work was followed up by Ashwell (1962), where 
localisation of stresses at two diagonally opposite corners is found for a strip in its 
first buckling mode. The actual shape of the strip was not computed. 

The constraint of zero Gaussian curvature causes the buckling patterns of inex- 
tensible strips to differ strongly from the oscillatory buckling patterns of extensible 
strips. The helicoidal shape of the edges of the strip found in Coman & Bassom 
(2008) is not a solution for inextensible strips not even for infinitesimally small 
twist. Rather we observe a sequence of relatively flat triangular domains which are 
not restricted to the core of the strip. The edges show a sequence of points with 
high curvature. 

It is worth noting that both responses, the smooth sinusoidal one and the lo- 
calised one, can be observed on the same paper strip model depending on the envi- 
ronmental conditions. When the humidity is low and the paper is dry it behaves like 
an inextensible material. When it is slightly wet it becomes noticeably stretchable. 
Note also that the solution for a slightly extensible strip is well described by the 
inextensible model almost everywhere except for small domains where the stress 
concentrates. In this paper we compute geometrically-exact developable solutions 
of inextensible strips. 

A geometrically-exact set of equilibrium equations for the large deformation 
of thin inextensible plates of finite width has recently been derived (Starostin & 
van der Heijden 2007). The equations are ordinary differential equations and ob- 
tained by using the inextensibility constraint to reduce stresses and strains to the 
centreline of the strip. They are much easier to analyse than the usual partial dif- 
ferential equations of elastic plate theory. The new set of equations was used to 
solve a classical problem in mechanics, namely to find the equilibrium shape of a 
developable Mobius strip (Starostin & van der Heijden 2007). Numerical solutions 
revealed the existence, for any aspect ratio of the strip, of a nearly flat triangular 
region associated with the (unique) inflection point of the strip (see figure [Ijb)). 
The triangular facet of the Mobius strip solution clearly resembles the facets of the 
buckling pattern of the twisted strip in figure [Ija) and in this paper we use the 
new system of equations to construct post-buckling solutions in good agreement 
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(a) (b) (c) 

Figure 1: (a) Twisted acetate model strip under tension, (b) Mobius developable 
structure of aspect ratio 27r with generators shown, (c) Trapezoidal segment taken 
from the Mobius structure to construct the periodic strip, with normal, binormal 
and centreline shown, ?7'(0) — 0,'q{L) = 0. The colouring changes according to 
the local bending energy density, from violet for regions of low bending to red for 
regions of high bending. Note the singularity on the edge of the strip. 

with experiment. The central idea is to modify the boundary-value problem for 
the Mobius strip such as to 'cut out' the triangular (more precisely, trapezoidal) 
region (see figure [Ijc)) and to use symmetry to reflect and multiply the elementary 
triangular facet into a periodic triangular pattern. This procedure avoids having to 
integrate numerically through the bending energy singularity associated with the 
vertex of the triangle on the edge of the strip, as found in Starostin & van der 
Heijden (2007) and analysed further in Hornung (2009). 

Having obtained a post-buckling strip solution with a certain number, say n, 
of triangular facets we then study the strip's force-extension and moment-twist 
behaviour for various mode numbers n. Neither gravity nor other distributed forces 
acting on the strip are taken into account in the present work. 

The paper is organised as follows. In Section 2 we formulate the boundary-value 
problem for the centreline-reduced equations for a developable strip. We also use 
symmetry properties of the solution to construct triangular buckling patterns for 
the strip by concatenating elementary facets. In Section 3 we numerically solve 
the boundary-value problem and compute load-displacement curves for the post- 
buckling strip solutions for various mode numbers n. In Section 4 we discuss our 
results and draw some conclusions. 

2. The boundary- value problem 

A sufficiently thin elastic surface will deform by bending only (Witten 2007) and 
therefore deform isomctrically. If such a surface is flat in its unstressed state it will 
remain so under deformation and therefore have zero Gaussian curvature (Graustein 
1966). It is said to be developable. 

Now consider a rectangular sheet or strip. If r(s) is a parametrisation of the 
centreline of the strip, s being arclength, then 

x(s,i)=r(s)+t[b(s)+77(s)t(s)], 
t{s) = j]{s)k{s), s e [0, 1], t e [—w, w] 

is a parametrisation of an embedded developable strip of length I and width 2w 
(Randrup & R0gen 1996). Here t and b are two unit vectors of the Frenet frame 



Article submitted to Royal Society 



A. P. Korte, E. L. Starostin and G. H. M. van der Heijden 



{t,n,b} of tangent, principal normal and binormal to the centreline, while k and 
T are, respectively, the curvature and torsion of the centreline, which uniquely 
specify (up to EucHdean motions) the centreline of the strip (Graustein 1966). By 
equation (I2.1|) . the surface, in turn, is completely determined by the centreline 
of the structure. The straight lines s = const, are the generators of the surface, 
which make an angle arctan(l/77) with the positive tangent direction. Given this 
parametrisation, the mean curvature M can be easily calculated, e.g., by using the 
coefficients of the first and second fundamental forms of the surface, themselves 
calculated from partial derivatives of x with respect to s and t (Graustein 1966). 
The result is 

M=-^l±^, (2.2) 

2 1 + tT]' ' ^ ' 

where the prime denotes differentiation with respect to arclength s. 

Introducing rectangular co-ordinates (ui,U2) by developing the surface into a 
rectangle, 

ui = s + trj{s), U2 = t, (2.3) 

the bending energy of a strip of thickness 2h can be written as the following integral 
over the surface of the strip (Love 1927): 

U = 2D Af 2 dui du2, (2.4) 

where D = 2Eh^/[3{l — v'^)] is the flexural rigidity, E is Young's modulus and 
1/ is Poisson's ratio. When M (given by (|2.2p ) is substituted into (j2.4p . and the 
co-ordinates changed to s, t, the t integration can be carried out (Wunderlich 1962) 
giving 

U = Dw f giK,r,,r,')ds, (2.5) 

Jo 



with 



,(«,,,,',. K»{l + fl'^>o!(l^). («) 



In the zero- width limit, w — > 0, this reduces to the Sadowsky result g = k^(1 + 77^)^ 
(Sadowsky 1931). 

Minimisation of this elastic energy functional is a one-dimensional variational 
problem cast in a form that is invariant under Euclidean motions. In Anderson 
(1989) Euler-Lagrange equations for such geometric variational problems are de- 
rived in a general context via a splitting of the cotangent bundle T*J°° of the 
infinite jet bundle J°° of a fibred manifold, which induces a bigrading of the dif- 
ferential forms on J°° known as the variational bicomplex (see also Kogan & Olver 
2003). In more physical terms they can be written in the form of six balance equa- 
tions for the (normalised) components of the internal force, F = {Ft,Fn, Ft,)'^ , and 
moment, M = {Mt, Mn, Mt,)'^ , in the directions of the Frenet frame, and two scalar 
equations (Starostin & van der Heijden 2007, Starostin & van der Heijden 2009): 

F' +UJX F = 0, M' + LJxM + txF = 0, (2.7) 

d^g + T]Mt + Mb ^ 0, (d^'g)' ~d^g~KMt^O, (2.8) 
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where lo = K{ri, 0, 1)^ is the Darboux vector. The equations (12.71) are nothing but 
the vectorial fixed- frame force and moment balance equations F' = 0, M' + r' x 
F = 0, written out in the Frenet frame. (Here we adopt the convention that bold 
roman symbols are used for vectors while bold italic symbols are used for triples of 
components of these vectors in the Frenet frame.) It follows immediately that F ■ F 
and M ■ F are first integrals of the equations. Note that the first equation in (|2.8p 
is algebraic in the variables (k, tj, tj'), while the second is a second-order ordinary 
differential equation (ODE) in rj. 

It will be of interest to a more general audience how the same equations can 
be obtained from first principles using standard variational methods, extending 
to a function of k, t, kI and r' the examples in Capovilla et al. (2002) (which 
covers functional dependence up to «;,t only). This is done in the appendix, where 
it is shown that it is straightforward to accommodate the additional functional 
dependence on k',t' without computing any new variations, by simply reusing the 
results already given in Capovilla et al. (2002) and Langer & Perline (1991). The 
additional terms generated by the dependence of g on k' and r' are easily managed. 
It is straightforward to extend this method to a functional involving any number 
of derivatives of k and r, and hence obtain the closed-form expressions given in 
Starostin & van der Heijden (2009). 

The shape of the strip's centreline is found by first differentiating the first equa- 
tion in (12.81) in order to turn it into a differential equation and then numerically 
solving equations (|2.7p and (|2.8p as a boundary-value problem (BVP) in conjunc- 
tion with three Euler- angle equations describing the evolution of the Frenet frame 
relative to a fixed frame, and the centreline equation r' = t. Taking Love's con- 
vention for the Euler angles (0, V', 0) (Love 1927), the derivatives of the angles are 
related to the curvature and torsion as 

& = Kcoscf), 

ip' = K sin (/)/ sine, (2.9) 

(j)' — —K sin (j) cot 9 + KT], 

while the centreline equation in component form gives 

(2.10) 



Equation (j2.9p can be written down directly from the Darboux vector by noting that 
Love's convention is the usual y-convention (van der Heijden & Thompson 2000) 
in classical mechanics with ip and interchanged (Goldstein 1980). Alternatively 
it can be obtained from the Frenet-Serret equations. The convention is such that 
the polar singularity (here at ^ = 0) usually associated with Euler angles does not 
cause any problems for the solutions we are interested in (the fiat strip will have 
9 — 7r/2). Alternatively a 4-parameter quaternion representation can be used to 
avoid the singularity at the expense of a norm condition. 

Before we specify the boundary conditions let's have a closer look at the buckling 
pattern in figure [TJa) and the Mobius strip solution in figure [TJb). The Mobius strip 
solution has special points where either rj or r/ is zero. Points where 77' = are called 
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cylindrical points as the surface is locally a cylinder, i.e., the mean curvature in 
is constant along the local generator. In figure [IJb) this corresponds to a generator 
of constant colour. Points where 77 = 0, by contrast, are called conical because the 
edge of regression, on which nearby generators intersect each other, has a cusp. At 
these conical points the generator is perpendicular to the centreline, as follows from 
(|2.ip . Clearly there must be at least one point where 77' = between any two points 
with 77 = 0. In fact, the Mobius strip has three cylindrical points and three conical 
points, one of the latter being special because it corresponds to the only inflection 
point of the centreline (where k = 0). At this point the binormal component of both 
the force and the moment are zero. Furthermore it is found that at the inflection 
point, \r]'\ — J> 1/w (i.e., the edge of regression reaches the edge of the strip) and 
the bending energy density g diverges, i.e., we have stress concentration. As is seen 
in figure [IJb), coming out of this singular point is a nearly fiat (violet) triangular 
region. 

Now, turning to figure [Ha) we observe that the buckling pattern consists of 
points of high stress located alternatingly on both edges of the strip, while locally 
cylindrical ridges bound flat triangular (more precisely, trapezoidal) regions similar 
to those found in the Mobius strip solution. This suggests that we can describe the 
buckling pattern by a solution of the equations built up of alternating copies of 
the trapezoidal section between the inflection point (where 77 = 0) and the nearest 
cylindrical point (where 77' = 0). A cut-out of this section is shown in flgure[TJc). 
Note that both bounding generators are of constant colour, illustrating that the 
section can be reflected about both end generators. 

A two-step symmetry operation is therefore used to construct a strip of length 
2nL from a trapezoid of length L. Let the arclength parameter of the centreline be 
s = at the cylindrical point (7/ — 0) and s = L at the inflection point {rj = 0) (see 
figure [TJc)). The first operation is a rotation through 180° of the trapezoid about 
the normal ni at s = (see also figure [5]). The original and the rotated trapezoid 
together make a continuous surface, which forms one period of the full strip (Pi of 
figure [5]), of length 2L. The binormal to the strip at the inflection point is denoted 
bo. 

In the second step a symmetrical strip of 7i -I- 1 periods is obtained from a parent 
strip of n periods by rotating an end period of the parent strip 180° about the end 
binormal. This binormal is located at the inflection point of that end period and 
therefore aligned with the end generator. Thus strips of any period can be built up 
in this way by successively reflecting an end period about its end binormal. The 
length of the centreline of a symmetrical strip of n periods made in this way is 
therefore 2nL. In particular, let r(0) = Tq, r{L) = ri be respectively the position 
vectors of the centreline at the inflection and cylindrical points of the trapezoid 
solution of the BVP. Deflne a rotation of tt around the unit vector g, 

i?g(a) = 2g(g.a)-a, (2.11) 

then the flrst and subsequent periods are obtained by the iteration 

r2 =ri +i?„i(ro -ri), b2 = i?„i(bo), 

r2i = r2i-2 + ^b2,-2('"2i-4 — r2i-2), b2i = -Rb2i-2(b2i-4), 



(2.12) 



where i — 2, . . . ,n. The h2i are then the unit binormals at the inflection points of 
the periods of the strip. The reflection rules in (j2.12p are for the centreline, but 
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Figure 2: Mode n = 4 synimetry operations. The first half period is rotated about 
til to form the first period Pi. Rotating Pi about b2 produces the second period 
P2. Rotating P2 about b4, then P3 about bg, produces the final two periods of the 
figure. 

can be easily generalised to the whole surface, as was done in Starostin & van der 
Heijden (2007) to obtain closed-strip solutions from the trapezoid solution of the 
BVP. These refiection rules are shown in figure [5] for the example n — 4 with the 
relevant binormals shown. The initial trapezoid (with one boundary aligned with 
bo) is rotated about rii giving the first period Pi. This first period is then rotated 
about b2 to give the second period P2 , which itself is rotated about b4 to give the 
third period P3. Finally the third period is rotated about bg to give the fourth 
period P4. 

It remains to formulate the boundary conditions for the initial trapezoid. At 
s — these are 

F„(0) = 0, Af„(0) = 0, v'iO) = 0, .2 13) 

x(0) = 0, y(0) = 0, z(0) =0, 



where the values for x, y and z fix an arbitrary position in space, and 
Afb(O) = -r^{0)Mt{Q) - 2k(0)(1 + ry2(0))2, 



(2.14) 



which can be read off as the s = limit of the first equation of (|2.8p . where ri'{Q) = 0. 
Equation (j2.14p is required to fix the integration constant for the ODE obtained 
by differentiating the first equation of (|2.8p . (Note that by Taylor expanding the 
second of equations (|2.8I) around s = 0, one can also show that M((0) = 1(1 + 
772(0))k(0)(-67?(0) + w;2(i + Tf{0))r^"{Q))). 
The boundary conditions aX s = L are 



k(L) = 0, Ffc(L) = 0, Mfc(L) = 0, 
0(L) = 7r/2, V(i) = 0, 0(L)=^, 



(2.15) 



where the angles fix an arbitrary orientation of the strip in space. 

The force and moment boundary conditions in (j2.13p and (|2.15p are enforced by 
the rotational symmetry. Since rii and bo are local axes of rotational (or refiection) 
symmetry, any component of force and moment along those axes must vanish. This 
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guarantees continuity of forces and moments in the strip and therefore yields a vahd 
solution of the mechanical problem. 

The remaining boundary conditions for the system of ODEs are the projections 
of the end force and moment of the strip on the end-to-end unit vector e = e/|e|, 
e = r2„ -ro, 

F(i)-e = F, M{L)-e^M, (2.16) 

where F, M are also continuation (control) parameters. Equations (j2.13D - (j2.16p 
give a final count of 15 boundary conditions for the same number of first-order 
ODEs. 

The BVP is solved by continuation of each of the boundary conditions in (|2.16p 
in turn. To obtain the starting solution of the continuation, the numerical solution 
of the Mobius strip between the inflection point and nearest cylindrical point is 
continued in i^„(0) and A/„(0) until equations (|2.13p are satisfied. 

There are significant numerical difficulties solving this BVP as both ends of the 
integration interval have a singularity. The boundary condition k{L) = enforces 
an inflection point at x ~ L. At this point the torsion r must also be zero and in 
fact it must go to zero, in s, faster than the curvature n (Randrup & R0gen 1996). 
Therefore we also have ri{L) = 0. In practice, to compute a starting solution we first 
compute an approximate solution with k{L) ~ 0.1 to stay away from the singularity 
at X = L. When all other boundary conditions are satisfied we 'pull' the solution 
into the singularity by continuing k{L) to zero as far as possible, typically reaching 
values of 0.001. At this point we typically have ri{L) ~ 0.002, while |f7'(i)| — l/w, 
the distance from the singularity, is typically as small as 10^^. 

At the other singularity, at s = 0, numerical convergence requires Taylor ex- 
pansions (up to fourth order in our case) of the left-hand sides of equations (12. 8p 
about ry' = to be used for a small interval around s = 0. In addition, to improve 
convergence, the absolute value of the logarithm is taken in equation (j2.6|) . We 
note that it has not proved possible (for instance by similar Taylor expansions) to 
remove the singularity in (j2.8p at rj' = 1/w. 

As a check on the numerical results, the first integrals F ■ F and M ■ F are 
typically found to be constant to within 10~^. 

Once the BVP has been solved, the surface of the strip is obtained from (12.11) 
and the refiection rules (|2.12[) . We note that no measures in the above formulation 
are taken to prevent self-intersection of the strip. 

3. Numerical results: response of the strip to applied loads 

The AUTO continuation software (Doedel et al. 2007) was used to compute response 
curves of applied force F = F(L) • e against end-to-end distance of the strip |e|, 
and applied moment M = M(L) • e against the accumulated twist angle a. As 
the component of bo perpendicular to e is ao = bo — (bo • e)e, with a2„ similarly 
defined, the twist angle a is the angle between the components of the two end-of- 
strip binormals perpendicular to the end-to-end vector of the strip, and therefore 
cosa = ao • a2„/(|ao||a2„|). 

Under a distance rescaling s —!■ s' = sk, the force and moment scale as F — ;> F/fc^ 
and M -^ M/fc (M scales as k, cf. equation (j2.14p ). k scales as k — )► K/k, whereas 
7] and the Euler angles are scale invariant. Rescaling can be used to obtain strip 
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(a) 




(b) 




Figure 3: Shapes of strips for n — 8 
M = 2.71. (c) F = 0.197, M = 12.15. 



(c) 
(a)F 



13.42, M = 0.162. (b) F = 6.80, 



solutions of a given period from solutions of another different period, with the same 
aspect ratio, via a w continuation. For example, starting from a strip of length 2nL, 
width 2w, and period n — 2m, selecting n/2 consecutive periods gives a strip of 
length nL, continuing this solution to w and rescaling s ^^ 2s gives a period n/2 
strip with the same aspect ratio as the original strip. All strips are scaled to an 
aspect ratio of 2nL/2w — 10.5300 unless otherwise stated. 

Many different strip shapes are obtained depending on the boundary conditions 
(j2.16p . In figure [3] arc shown three shapes of a strip with n — 8, one under a high 
axial tension, one under a relatively high axial moment and an intermediate one 
which is in reasonably good agreement with the experiment in figure [IJa). Here, 
and in all surface plots in this paper, the colouring varies according to the bending 
energy density, from violet for regions of low bending to red for regions of high 
bending (scales are individually adjusted). 

Continuation results for modes n = 2, 4, 8 are shown in figures [4H6] (for aspect 
ratio 2nL/2w — 10.5300, except for figurelHUa), which has aspect ratio 42.30). On 
the left of each figure is shown the force response for the sequence of axial moments 
(M ) shown on the legend. The length scale is arbitrary and is such that L = 0.6581. 
On the right of each figure is shown the moment response for the sequence of axial 
forces (F) shown on the legend. These curves were obtained by continuation, where 
F was varied, keeping M fixed or vice versa. Exceptions to this are some curves in 
figure [ni[a,b) and figure [Sfa), where both F and M were allowed to vary, only so 
that a continuation could be started at an arbitrary point on the graph (at which 
point one of F, M would then be held fixed). This accounts for the curve-crossing 
in this figure. 
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2 4 6 8 10 

e 

(a) (b) 

Figure 4: n = 2 mode, (a) Force F versus scaled end-to-end distance, jej. The curves 
of constant M form a nested sequence on the plot, with M increasing from the top 
middle dashed curve (M = —4.15) to the dotted curve, then from the dot-dashed 
curve to the dashed curve on the right, (b) Moment M versus twist angle a (in 
radians). The curves of constant F form a nested sequence on the plot, with F 
increasing from the dashed curve {F — —0.9119) to the dotted curve. 





(a) (b) 

Figure 5: n — A mode, (a) Force F versus scaled end-to-end distance, |e|. The curves 
of constant M form a sequence on the plot, with M increasing from the top left 
dashed curve (M — —25.1) to the dotted curve, then from the dot-dashed curve 
to the top right dashed curve, (b) Moment M versus twist angle a (in radians). 
The curves of constant F form a sequence on the plot, with F increasing from the 
dashed curve {F = —4.06) to the dotted curve {F = —3.26), then to the top dotted 
curve. 
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Figure 6: n = 8 mode, (a) Force F versus scaled end-to-end distance, |e|. The 
curves of constant M form a nested sequence on the plot, with M increasing from 
the left dashed loop {M — —1.34) to the dotted curve, then from the dot-dashed 
curve to the top dashed curve, (b) Moment M versus twist angle a (in radians). 
The curves of constant F form a nested sequence on the plot, with F increasing 
from the dashed curve {F = —12.87) to the dot-dashed curve. 



't^ 150 




-50 



M=-9.19 
M =13.67 
M=-35.13 







Figure 7: n — 8 mode. Force F versus scaled end-to-end distance, |e|, showing strip 
solutions. 



In figure m^a) , for the n = 2 mode, the sequence of curves at constant negative 
applied end moment is the group at the bottom left, with the moment increasing 
from top (large negative) to bottom (near zero) . The group of curves at the bottom 
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(a) (b) 

Figure 8: Shapes of strips, (a) n — 2, figure-of-eight, (h) n = 8, F — —19.975, 
M = 0.1424. 



right of the figure is for constant positive applied end moment, with the moment 
increasing from bottom (near zero) to top (large positive). In figure Ufb) is a se- 
quence of curves at constant applied end force, increasing from bottom (negative) 
to top (positive). The curves for smallest (negative) constant force start at zero 
moment, have a maximum twisting angle, which returns to zero twist angle, which 
is the figure-of-eight shown in figure [SKa). 

In figure [HJa), for the n = 4 mode, the sequence of curves at constant negative 
applied end moment is the group at the bottom left, with the moment increasing 
from top (large negative) to bottom (near zero). The group of curves at the bottom 
right of the figure is for constant positive applied end moment, with the moment 
increasing from bottom (near zero) to top (large positive). In figure [5Kb) is a se- 
quence of curves at constant applied end force, increasing from bottom (negative) 
to top (positive) , showing the response curve of applied end moment to the twist 
angle a. The three lowest negative applied force continuations terminate at a = 27r, 
again forming a closed strip, this time the T^^i elastic torus ribbon knot. 

In figureEKa), for the n = 8 mode, is a sequence of curves at constant applied end 
moment (positive and negative). This picture is more complicated than the other 
modes. Some of the curves terminate on the vertical axis, i.e., when the end-to-end 
distance vanishes. These closed ribbon solutions are invariably torus ribbon knots. 
For example, the upper dashed curve {M = —20.816) terminates in the double 
cover of the T^^i elastic torus ribbon knot. The series of nested closed curves have a 
quadruple cover of the figure-of-eight solution where the curves approach vanishing 
end-to-end distance. In figure EKb) is a sequence of curves at constant applied end 
force, increasing from bottom (negative) to top (positive), showing the response 
curve of applied end moment to the twist angle a. 

In figure [71 for the n — 8 mode, is a sequence of curves at constant applied 
end moment, with some of the resulting structures shown. The upper curve (M = 
—35.13) terminates in the Tg_i elastic torus ribbon knot. Both ends of the dot- 
dashed curve (M = 13.67) terminate in a double covering of the T4^i elastic torus 
ribbon knot, whereas the Tg 3 structure (shown at the minimum end-to-end distance 
of the loop) does not quite close, as can be seen by the fact that the curve from 
which it originates does not meet the vertical axis. The other structure shown on 
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Figure 9: Fold structure of scaled end-to-end distance versus applied force F, n — 
2, . . . , 8 at constant M = 14.49. 



the same loop in the line diagram (shown at the maximum end-to-end distance of 
the loop on the figure) is related to the Tg^a torus knot, but it is also not closed, 
and moreover the surface of the strip is self-intersecting. 

As an example of further strip solutions, in figure [HIb) is shown a solution, 
with n — 8 and under a compressive axial force, that is not located on any of the 
computed response curves. 

Figure [S] displays force-extension curves for varying mode number, n = 2, . . . , 8, 
at fixed moment M. To obtain this plot, strips for different modes have been scaled 
to the same aspect ratio 2nL/2w = 10.5300. Parts of these curves with positive 
slope are expected to correspond to stable solutions. The curves predict that, for 
n > 5, under increasing tension solutions jump to higher mode (down in |e|) as the 
force is increased beyond the folds seen in the diagram. 



4. Discussion 

When twisted, and pulled, an acetate model strip buckles into a regular pattern of 
triangular facets. We have computed periodic solutions describing this buckling pat- 
tern by formulating and numerically solving a geometrically-exact boundary- value 
problem for the large deformation of a thin, wide, inextensible strip. We have also 
obtained response curves of force against end-to-end distance and twisting moment 
against end-to-end angle for mode numbers n = 2, 4 and 8. Our results predict 
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critical forces and jumps into higher buckling modes that would be interesting to 
explore experimentally. 

By construction our solutions are periodic, which tends to be what one observes 
in experiments. However, non-periodic solutions can be constructed in a similar 
way by matching different trapezoidal segments. One would keep the first symme- 
try operation (reflection about the normal at the cylindrical point) but instead of 
reflecting about the binormal at the inflection point in the second step one could 
match the segment to a (suitably rescaled) trapezoidal segment of different length 
L. The resulting solution would no longer be symmetric about the inflection points. 

In Mansfield (1989) and Ashwell (1962) analyses were performed similar to ours 
in that there too the double strain energy integral was reduced by integrating along 
the generator. However, the integrand expression was subsequently simplified by 
using an approximate moment balance equation that ignores the non-planarity of 
the strip. This approach gives a realistic evolution of the generator with arclength, 
but fails to predict the appearance of a cylindrical point at finite distance from the 
clamped end of the strip. 

Triangular patterns are known to occur in a variety of problems of elastic sheets, 
including fabric draping and paper crumpling (Witten 2007). A sheet crumples 
(forming sharp points and straight creases) when it is forced into a constrained 
area. The sheet is predominantly under compression. It is interesting to note that, 
by constrast, in our triangular buckling pattern the strip is in (relatively high) 
tension. In both cases we observe a focussing of the strain energy which may lead 
to fracture of the material. Strain energy localisation thus appears to be a generic 
response of a thin elastic sheet to an external constraint. 

Our results may be relevant for paper, fabric and sheet metal processing. They 
may also be of importance in the robotic manipulation of flexible belts, e.g., film 
circuit boards (Wakamatsu et al. 2007). In all these sheet manipulations it is im- 
portant to avoid shapes with high concentrated stresses that may lead to tearing. 
Our formulation may help to choose boundary conditions that avoid unwanted 
configurations. 

This work was supported by the UK's Engineering and Physical Sciences Research Council 
(EPSRC) under grant no. EP/F023383/1. 



Appendix A. Derivation of the equations using standard 

variational techniques 

We wish to calculate the variation of the functional H — / /(k, r, k',t') ds. This 
calculation can be done using the formulation of Capovilla et al. (2002), substituting 
explicit expressions for quantities wherever they arise and grouping similar terms, 
which gives fully simplified expressions for all quantities of interest. 

Alternatively, one can re-use the variations already calculated in Capovilla et 
al. (2002), to yield the same results, without having to compute new variations. 
There, for example, the variation is given for Hi = J /i(k, r) ds, stating that it can 
be written down from the variations of H2 — J /2(k) ds and H3 — J fsi^r) ds, which 
are previously calculated. The resulting force and moment corresponding to Hi can 
then be written down from the expressions for the force and moment for H2 and 
H3. This follows from the chain rule. Similarly, one can calculate the variation of 
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the functional H — J /(k, r, k', r') ds using SHi without computing new variations. 
This follows as a consequence of the chain and product rules in differentiation. The 
resulting force and moment corresponding to H can then be written down from the 
expressions for the force and moment for Hi, with a small number of additional 
terms arising from the chain and product rules. 

Instead of following these approaches, we here use a third and more economical 
method by leaving SH expressed in terms of Sk and St only. This method is easily 
generalised to functionals involving higher derivatives of k and r. Following the 
notation of Capovilla et al. (2002), the infinitesimal deformation of a space curve 
is Sr = ^||t + \l/in + ^2b, and denoting 5\\,d± as the tangential and normal parts 
of the deformation respectively, the variation of any functional H = J f ds is 

SH ^5aH+ f 5j_fds, (Al) 

where 6aH = 5\\H + Jf5j_ds = J ((/*||)' - /k*i) ds. 

We therefore need to calculate the second term on the right-hand side of equation 
(|A 1[) . where by the chain rule, 5±f = f^S^K + frS^r + /k/S^k' + fr'S^r' , using 
the notation /^ = df/d/3. Since for any scalar h, d±{h') = k/i'^i + {d±hy , then 

S±f - {f.-fL')S±'i+{fr-K,)S±T + K{K'U,+T'fr')^l + {US^K + fr'S^Ty, (A 2) 

where 5±k and S±t are given in Capovilla et al. (2002). Equation (|A 1[) along with 
equation (jA2p is of the form 



SH= / ds£,*, 



JdsQ', (A3) 



with Q the Noether charge, so that the Euler-Lagrange equations £i — can be 
immediately written down, since these are just the coefficients of ^^ which can 
be read off from equation (jA ip . These are two coupled differential equations in 
the unknowns k,t, and the known density /, which are quoted in general form in 
Thamwattana et al. (2008) and in Hangan (2005) for the Sadowsky functional. 

Apart from the Euler-Lagrange equations, it is of interest to find expressions for 
the conserved force F and the moment M. The force F is obtained by specialising 
the deformation to a constant infinitesimal translation 6r = e, where F is defined 
by Q = — e ■ F. Similarly, the conserved moment T is obtained by specialising the 
deformation to a constant infinitesimal rotation Sr = ft x r, where T is defined by 
Q = — J7 • T and decomposed into T = r x F + M. Thus only the total derivative 
parts of (jA ip contribute to F and M. But these contributions can be read off from 
the results in Capovilla et al. (2002) by noting that a term in equation (|A2p of the 
form aS±b, where b = k,t, has to be re-expressed in the form (jA 3p to isolate its 
total derivative part. This has already been done in Capovilla et al. (2002). For a 
term in equation (jA 2p of the form {a5±b)' , one just needs the contribution from 
Si_b. 

Thus, for an infinitesimal rotation (5r' = il x r' = fi x t and 5t" — kJI x n. 
Using the compact expressions of Langer & Perline (1991) (equations (3.7b,c)), 5k = 
Sr" ■n—2KSr' -t gives Sn = 0, i.e., (5_lk = —5\\k = — ^||«;'. But ^|| = Sr-t = (rxt)-J7, 
giving (5_lk = — J7 • {n'r x t). As this is of the form — n ■ (r x F), Q = aSj^n in (|A2p 
does not contribute to the moment M. Similarly, St — Sr"-h/K'+5r'-(Kh — Tt) = 0, 
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giving S±T — —ft ■ (r x r't). As this is of the form — $7 • (r x F), Q = aS±T in (jA2p 
also does not contribute to the moment M. (Note that by retaining these terms one 
would obtain the force F, but in the next paragraph we will, instead, calculate F by 
considering a constant infinitesimal translation.) Now, from Capovilla et al. (2002) 
(equations (49) and (67)), a6±K contributes —ah to the moment, whereas a6±T 
contributes — at — a'n/K. Thus from equations (|A II) and (|A2[) one can immediately 
write down the moment as: 



M = (/;, - fr)t + if:!, - f^)n/n + if,, - /Jb. (A 4) 

In order to calculate the contribution to the force of Q = fK'S±K + Jt'S^t, one 
notes that since (5r = e is a constant, it follows from the formulation of Langer 
& Perline (1991) (equations (3.7b,c)) that 5k — 0, i.e., 5±k = —S^^k = — ^||k' = 
— e • (/i't). Thus the contribution of Q = f^'^^n in (IA2p to F is K'/^'t. In a similar 
manner, for this constant translation, 5t — 0, giving 5±^t — e ■ (r't), so that 
Q = Jt'S^t in (jA2l) contributes r'/^'t to the force. Now from Capovilla et al. 
(2002) (equations (48) and (64)), a5±Hi contributes nat + a'n + rah to the force, 
whereas a5±T contributes rat + tu'vi/k — {{a' / k)' + Ka)b. Also, 5\\H contributes 
— /t to the force. Thus from equations (|A ip and (|A2I) one can immediately write 
down the force as: 

F^{-f + k'U + T'fr' + K{f, -/:,)+ <fr - f'r')) t 

+ (/^-/:, + ^(/;-/;))n (A5) 

+ (t(/k - /:0 - ^{fr - fr-) - ((/; - /"0/^)')b. 

From these components of F and M , the differential equations for the moment 
components and the differential equation for Ft in equation (|2.7p are easily obtained. 
For instance, from (jA4[) it is seen that M[ — kM„, while equations for M„, Mf, 
and Ft can be extracted similarly. The remaining two equations for Fn and F^ are 
obtained from (jA5|) and the Euler-Lagrange equations 8i = Q. 

In summary, the variation oi H = J /(k, t, k', t') ds can be easily obtained from 
the variations Sk,St. Their explicit expressions are not required; their contributions 
can be read off from Capovilla et al (2002). The additional terms generated by the 
additional dependence of / on k' and r' are therefore easily managed. 

One way of regarding equations (|A4I) and (IA5|) . is that they prescribe F and M 
once K and r for the given density / are known. The two Euler-Lagrange equations 
for K and t are given by £i = 0, i.e., by setting the coefficients of ^^ to zero in 
equation (jA 3[) . Instead of solving the problem this way, however, it is preferable 
to set up an equivalent system of coupled one-dimensional ODEs as in Starostin & 
van der Heijden (2009). One way to do this is to use the natural variables F and 
M, and, for the functional p.6p . to change variables from {k,t) to {K,ri). 

By considering functionals of the form g{K,r],r]') = f{K,T,K',T'), one can show 
that {dg/dK)r,.r,' = /k + vfr + v'fr', {dg/dri),^r^> = nfr 4- k'/t' and {dg/dr]'),^^, = 
Kfr'- Using these identities, and writing the components of the internal force F 
and moment M in the directions of the Frenet frame of tangent, principal normal 
and binormal as F' = {Ft, Fn,Fi,)'^ , M — {Mt, Mn, Mi,)'^ , one can show, using the 
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tangent and binormal components of M from equation (jA4p , that 

d,g + r]Mt + Mb = 0, ^^ ^ 

, (A 6) 

{d,fg) ~ d,jg ~ nMt == 0. 

The two scalar equations, (|A6[) . along with the six balance equations (equation 
(|2.7p ) for the components of the internal force F and moment M (see Starostin & 
van der Heijden 2007, Starostin & van der Heijden 2009), are then the equations 
satisfied by the centreline of the developable strip for the unknowns F, M, k, tj: 

F' + UJ X F = 0, M' + LJxM + txF = 0, (A 7) 

d^g + ijMt + Mb - 0, [d^^g)' - d^g - nMt = 0, (A 8) 

where u) — ^(77,0, 1)"^ is the Darboux vector. Equations (jA7l) . (IA8|) constitute a 
system of differential-algebraic equations that are turned into a system of ODEs by 
differentiation of the algebraic equation in 
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